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An  Analytical  Approach  to  the  Selection  of  the  Reaction  Coordinate 
in  the  Calculation  of  Reaction  Profiles 

ABSTRACT 

An  algebraic  method  for  the  analytical  determination  of  an  expression 
for  the  reaction  coordinate  in  an  arbitrary  potential  energy  surface  is 
presented.  The  method  is  applied  to  the  thionyl  imide  -  thiazyl  S  hydrozide 
isomerization  reaction.  The  extent  to  which  the  method  leads  only  to 
required  reaction  coordinate  is  dependent  upon  the  manner  in  which  the 
potential  energy  surface  is  sampled  to  obtain  input  data  for  the  calculations. 
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INTRODUCTION 


Present  semi  empirical  molecular  orbital  methods  based  upon  the  modified 

Neglect  of  Differential  Overlap,  i.e.  MNDO  and  MNDO/3,  enable  the  geometry  of 

a  nwlecule  in  its  electronic  ground  state  to  be  accurately  calculated. ^  The 

method  also  provides  a  reasonable  accurate  estimate  of  the  standard  enthalpy 
(21 

of  formation.'  '  As  such  the  MNDO  method  can  be  utilized  for  the  calculation 
of  the  enthalpy  cf  activation  along  a  reaction  path  and  thus  provides  thermo¬ 
dynamic  information  about  a  reaction  mechanism.  This  is  generally  accomplished 
by  the  selection  of  a  "reaction  coordinate"  and  performing  a  series  of 
calculations  for  various  values  of  this  "reaction  coordinate"  being  certain 

that  the  calculation  connects  the  appropriate  reactant  and  product  states.  This 

(31 

is  shown  symbolically  in  the  fugure  below:'  ' 
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1.  One  generally  has  no  a  priori  knowledge  of  the  geometry  of  the  transition 
state  and  must  thus  make  no  assumption  concerning  its  geometrical  structure. 

Each  assumption  about  geometry  limits  the  class  of  reaction  path  to  be 
examined. 

(5) 

2.  One  frequently  encounters  energy  maxima  which  are  not  true  transition  states.'  ' 

3.  A  large  number  of  variables  have  to  be  examined  in  a  systematic  manner. 

This  generally  implies  that  a  large  number  of  calculations  must  be  carried 
out.  The  procedure  can  be  systemized  somewhat  by  the  use  of  "grid 
searches, but  the  number  of  calculations  required  is  still  far  in  excess 
of  what  would  be  required  if  one  had  direct  knowledge  of  the  reaction 
coordinate. 

The  principal  difficulty  in  this  type  of  calculation  is  the  recognition 
of  what  to  use  fo'"  the  reaction  coordinate.  The  energy  calculations  are 
performed  in  terms  of  a  set  of  cartesian  type  coordinates,  which  are  determined 
from  a  set  of  internal  coordinates  in  turn  composed  of  the  bond  lengths,  the 
bond  angles  and  the  dihedral  angles  necessary  to  specify  the  geometry  of  the 
molecule(s).  Since  the  geometry  of  the  transition  state  is  not  known  a  priori, 
the  transformation  which  relates  the  reaction  coordinate  to  the  internal 
coordinates  is  not  known. 

Furthermore  experience  has  shown  that  the  reaction  coordinate  does  not 
correspond  to  a  single  internal  coordinate  or  a  linear  combination  of  two  or 
three  internal  coordinates  but  is  rather  more  complicated. At  the  present 
time  reaction  coordinates  are  identified  by  a  "trial  and  error"  procedure, 
wherein  calculations  are  performed  by  variation  of  one  or  more  coordinates 
in  the  internal  basis  set  until  the  transition  state  is  found.  The  procedure 
is  very  inefficient  because  of  the  existence  of  other  maximum  on  the  multi¬ 
dimensional  energy  surface,  and  the  trial  and  error  method  in  fact  becomes 
intractable  for  systems  which  require  a  large  number  of  coordinates.  The 
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development  of  a  method  which  would  lead  to  a  systematic  identification  of  the 
reaction  coordinate  would  considerably  enhance  our  ability  to  perforin 
calculations  of  reaction  profiles. 


f 


As  pointed  out  in  the  introduction,  one  of  the  principal  difficulties 
encountered  in  the  calculation  of  the  profile  of  a  chemical  reaction  is  the 
recognition  of  what  to  use  for  the  reaction  coordinate  on  a  given  reaction 
mechanism.  The  reaction  coordinate  is  generally  not  a  simple  atom-atom 
distance,  bond  angle  or  dihedral  angle  but  it  is  expressible  as  a  linear 
combination  of  these  atom-atom  distances,  bond  angles  and  dihedral  angles. 
Furthermore,  it  can  be  expected  to  contain  those  coordinates  which  change 
markedly  as  one  proceeds  from  the  reactant  state  to  product  state.  The  actual 
reaction  coordinate  can  be  identified  and  an  analytical  expression  obtained 
for  it  by  the  method  outlined  below. 

A  comparison  of  the  geometrical  structures  of  the  reactant  and  product 
molecules  enables  one  to  identify  which  internal  coordinate  undergoes  a 
marked  change  as  one  passes  from  the  reactant  state  to  the  product  state. 

The  reactant  and  product  molecules  must  of  course  be  defined  using  the  same 
set  of  internal  coordinates  (i.e.  bond  distances,  bond  angles  and  dihedral 
angles).  This  can  always  be  done.  Let  the  subset  of  internal  coordinates  so 
identified  be  denoted  by  jq^i-  Suppose  n  such  coordinates  are  required.  The 
reaction  coordinate,  can  be  expressed  as  a  linear  combination  of  the  lq.jl, 
i.e. 

n 

0  =  f  a.q.  (1) 

Explicit  determination  of  the  reaction  coordinate  requires  a  knowledge  of  the 
coefficients  [a^j.  To  a  good  approximation  the  energy  of  the  system  is  a 
parabolic  function  of  the  reaction  coordinate  and  we  can  write 

E(())  =  A^^  +  +  C  (2) 

This  is  certainly  true  in  the  vicinity  of  the  transition  state  and  probably 
valid  in  regions  not  too  far  removed  from  the  transition  state.  Since  we 


I 


are  always  seeking  a  transition  state,  E  is  a  maximum  for  such  a  state  and  the 
calculus  tells  us  that  A<0  and  B>0. 

An  arbitrary  selection  of  k  sets  of  values  for  the  jq-l  enables  one  to 
calculate  k  values  of  the  energy,  E.  If  we  assume  the  energy  is  a  quadratic 
function  of  the  internal  coordinates,  [q|,  we  can  write 


E  =  ■£.  Z  a  .q.q.  +  Z  b.q  +  C 

i  j  ij  I  J  j  J 


(3) 


Let  q-,  denote  the  value  of  the  internal  coordinate  q,  in  the  set  m.  The 
^Im  1 

above  calculations  enable  us  to  write  down  a  set  of  n  x  (n+3)/2  equations  which 
can  be  viewed  as  a  set  of  linear  algebraic  simultaneous  equations  in  the 
coefficients  la--|,  |b.|  and  c.  Namely, 

I J  J 


'  1  j  *  :  ^j^jk  *  " 

k  =  1,2,3  .  n(n+3)/2 


(4) 


These  equations  can  be  solved  to  obtain  values  of  the  coefficients  la.  J, 

^  *3 

jb-l  and  c.  These  coefficients  are  in  turn  related  to  the  coefficients  of 
J 

equation  (1),  |aj.  Namely  using  equations  (1),  (2)  and  (4)  from  the 
quadratic  term  we  obtain: 


a. j  =  Aa.aj  (i.j  =  1,2  ...  n)  (5) 

from  the  linear  terms  we  get: 

b.  =  Ba.  (i  =  1 ,2,3  . . .  n)  (6) 

and  finally; 

c  =  C  (7) 


From  eqn.  (6)  we  in  effect  know  the  coefficients  |a^.  |  to  within  the  normalization 


constant  B.  If  we  redefine  the  norm  of  the  space  spanned  by  the  reaction 
coordinate,  we  can  calculate  the  jct-l  and  explicitly  obtain  our  reaction 
coordinate  from  Eqn.  (1). 

There  are  several  aspects  of  the  above  method  that  require  examination. 
They  are: 

(1)  The  question  of  the  effects  of  overcompleteness  and  undercompleteness 
with  respect  to  which  internal  coordinates  are  to  be  included  in  eqn. 

(1).  Does  undercompleteness  imply  failure  of  the  method?  Does  inclusion 
of  unnecessary  coordinates  cause  undue  algebraic  labor? 

(2)  The  extent  to  which  the  quadratic  assumption  made  in  eqn.  (3)  is  valid 
in  regions  of  the  space  far  removed  from  the  vicinity  of  the  transition 
state.  Another  functional  form  could  also  be  assumed,  i.e.  a  cubic  form 
for  E,  and  a  method  developed  from  it. 

(3)  The  use  of  various  norms  to  solve  eqn.  (6).  One  might  use  a  norm  based 
upon  0<^^_<  1 ;  a  norm  based  upon  the  range  of  distances  observed  in  the 
reactant  state  and  product  state  or  some  other  norm. 


RESULTS 


The  method  was  first  applied  to  the  isomerization  reaction  of  cis  thionyl- 
imide,  HNSO  to  cis  thiazyl-S  hydroxide,  HOSN.  This  system  was  selected  because 
this  reaction  has  been  studied  in  our  laboratory  and  the  potential  energy 
surface  is  well  understood  by  us,  including  the  detailed  structure  of  ttie 
transition  state  and  the  appropriate  reaction  coordinate(^^  The  structures  of 
the  reactant  and  product  are  shown  in  figure  1,  that  of  the  transition  state 
is  shown  in  figure  2.  The  reaction  coordinate  is  a  linear  combination  of  the 
N-H,  the  N-S  and  the  0-S  bond  distances.  The  internal  angles  do  not  play  a 
major  role. 

The  first  point  to  be  noted  is  that  the  method  becomes  very  inefficient 
if  the  number  of  internal  coordinates  entering  equation  1  becomes  very  large. 

In  the  table  below  are  listed  the  number  of  data  sets  required  (this  is  the 
same  as  the  order  of  the  set  of  linear  algebraic  equations  to  be  solved)  for 
a  specific  number  of  internal  coordinates  in  the  reaction  coordinate. 


Table 

Number  of  Internal  Coordinates  Which 
Enter  Reaction  Coordinate 

1 

2 

3 

4 

5 

6 

7 

8 
9 

15 

21 


Number  of  Data  Sets 
Required 

2 

5 

9 

14 

20 

27 

35 

44 

54 

135 

252 


It  should  be  recalled  that  a  fully  optimized  molecular  orbital  calculation 
must  be  performed  to  obtain  a  single  data  set. 

The  question  of  overcompleteness  and  undercompleteness  with  respect  to 
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which  internal  coordinates  must  be  included  in  equation  1  is  critical. 

Omission  of  a  critical  coordinate  will  lead  to  a  false  maxima  or  an  uttdersi rable 
transition  state.  The  inclusion  of  excess  internal  coordinates  in  equation  1 
will  cause  unnecessary  labor  in  the  calculation. 

The  results  obtained  for  the  reaction  coordinate  seem  to  depend  upon  tlie 
choice  of  data  sets.  If  the  data  sets  consisted  of  a  sampling  of  the  potential 
energy  surface  near  the  transition  state  the  method  yielded  the  correct 
reaction  coordinate.  If  the  data  sets  selected  spanned  a  wider  area  of  the 
potential  energy  surface,  the  method  generally  did  not  cor.  erge  directly  to 
the  reaction  coordinate.  This  appears  to  be  a  severe  limitation  on  the 
method.  This  may  be  due  to  the  limitation  of  the  use  of  a  quadratic  expression 
for  the  energy  in  terms  of  the  reaction  coordinate.  It  would  be  worthwhile 
to  investigate  other  functional  forms  for  equation  2. 

Since  the  project  only  lasted  for  a  brief  period  of  eight  weeks,  the 
number  of  points  that  could  be  investigated  was  limited.  The  question  of 
what  the  norm  to  use  for  the  "normalization"  of  the  reaction  coordinate  for 
example  was  only  superficially  examined.  A  norm  based  upon  ia.j|<l  for  all  i 
was  tried.  It  seemed  to  work  well  leading  to 

I)  =  .0897583  dQ_2  +  1.0  d5_^  +  0.186049  d^_^ 
for  the  reaction  coordinate. 

In  summary,  the  method  is  capable  of  yielding  the  correct  reaction 
coordinate.  It  must  be  judicially  applied  requiring  some  considerable  experience 
on  the  part  of  the  user  in  the  selection  of  the  data  sets  to  be  employed.  It 
possesses  a  real  value  in  that  it  can  be  used  to  rapidly  "zero  in"  on  a 
transition  state  of  the  user  has  some  a  priori  information  as  to  the  location 
of  the  transition  state  on  the  potential  energy  surface.  It  will  be  necessary 
to  examine  further  test  cases  to  ascertain  its  general  applicability  and 
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obtain  an  explicity  answer  as  to  whether  it  offers  an  efficient  alternative 
to  the  present  day  trial  and  error  procedure. 
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